##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(ggrepel)
library(patchwork)
library("scales")

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--mut_type"), type = "character"),
    make_option(c("--mut_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"

    sample_list_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"

    rsem_file <- "~/20220915_gastric_multiple/dna_combinePublic/mRNA/CombineTMM.DNAUse.NJMU_TCGA.tsv"
    mut_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/expression"
    gene <- "TP53"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
rsem_file <- opt$rsem_file
mut_file <- opt$mut_file
mut_type <- opt$mut_type
gene <- opt$gene

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_expression <- data.frame(fread(rsem_file))
dat_maf_public <- data.frame(fread( mut_file ))

##########################################################################################

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

dat_expression <- subset( dat_expression , gene_id == gene )
dat_maf_public <- subset( dat_maf_public , Hugo_Symbol == gene & Variant_Classification %in% Variant_Types )
colnames(dat_expression) <- gsub( "[.]" , "-" , colnames(dat_expression) )

## 分变异类型
# 错义突变
mis_mut <- c("Missense_Mutation")
# lof突变
lof_mut <- c("Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del" , "Splice_Site")

if(mut_type == "missense"){
    dat_maf_public <- subset( dat_maf_public , Variant_Classification %in% mis_mut )
}else if(mut_type == "lof"){
    dat_maf_public <- subset( dat_maf_public , Variant_Classification %in% lof_mut )
}

##########################################################################################
plotTpm <- function(tmm_combine = tmm_combine , out_name = out_name , title = title , width = width , y_tmm = y_tmm){

    if( length(unique(tmm_combine$Type_new)) == 3){
        my_comparisons_1 <- list( c(1, 2) , c(1,3) , c(2,3) )
        tmm_combine$Type_new <- factor( tmm_combine$Type_new , levels = c( "Mut\nShare" , "Mut\nNon-Share" , "Wild") , order = T )
    }else{
        my_comparisons_1 <- list( c(1, 2))
        tmm_combine$Type_new <- factor( tmm_combine$Type_new , levels = c( "Share" , "Other") , order = T )
    }
   
    y_max <- y_tmm
    if(y_max > 2000){

        if(y_max > 10000){
            y_breaks <- 10000
        }else if(y_max > 5000){
            y_breaks <- 5000
        }else{
            y_breaks <- 2000
        }
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.1) , breaks = c( 0 , 500 , 1000 , 2000 , y_breaks ) , trans = sqrt_trans() )
    }else{
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.1))
    }

    col <- c(
        rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
        rgb(red=2,green=100,blue=190,alpha=255,max=255) )

    plot <- ggplot( tmm_combine , aes( x = Type_new , y = TMM ) ) +
        geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6 , aes(color = Type_new)) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 1 , aes(color = Type) ) +
        scale_color_npg() +
        scale_fill_npg()+
        facet_grid(.~Class)+
        xlab(NULL) +
        ylab("TMM")+
        theme_bw() +
        ggtitle(title) +
        y_lab_lim +
        stat_compare_means(comparisons = my_comparisons_1) +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='right',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold',hjust=0.5,vjust=0.5),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.ticks.length = unit(0.2, "cm") ,
            axis.line = element_line(size = 0.5)) 
    return(plot)
}

plotTpm_combine <- function(tmm_combine = tmm_combine , out_name = out_name , title = title , width = width , y_tmm = y_tmm){

    trans <- function(num){
        up <- floor(log10(num))
        down <- round(num*10^(-up),2)
        text <- paste("p == ",down," %*% 10","^",up)
        return(text)
    }

    my_comparisons_1 <- list( c(1, 2) , c(1,3) , c(2,3) )

    result <- c()
    for( class in unique(unique(tmm_combine$Class)) ){
        tmp <- subset( tmm_combine , Class == class )

        a <- tmp[tmp$Type_new==unique(tmp$Type_new)[1],"TMM"]
        b <- tmp[tmp$Type_new==unique(tmp$Type_new)[2],"TMM"]

        p <- wilcox.test( a , b )$p.value
        if( p < 0.01 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "p == " , round(as.numeric(p) , 2) ) 
        }

        tmp$p_text <- ""
        tmp$p_text[1] <- p_text
        result <- rbind( result , tmp )
    }

    y_max <- y_tmm
    if(y_max > 2000){

        if(y_max > 10000){
            y_breaks <- 10000
        }else if(y_max > 5000){
            y_breaks <- 5000
        }else{
            y_breaks <- 800
        }
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.2) , breaks = c( 0 , 500 , 1000 , 2000 , y_breaks ) , trans = sqrt_trans() )
    }else{
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.2))
    }

    col <- c(
        rgb(red=0,green=159,blue=134,alpha=255,max=255) ,
        rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
        rgb(red=2,green=100,blue=190,alpha=255,max=255) 
        )

    result$Type_new <- factor( result$Type_new , levels = c( "Trunk" , "Private" , "Wild") , order = T )
    result2 <- result
    result2$Class <- "IGC+DGC"
    result <- rbind(result , result2)
    result$Class <- factor( result$Class , levels = c( "IGC+DGC" , "IGC" , "DGC") , order = T )
    
    print(unique(result$Type))
    plot <- ggplot( result , aes( x = Type_new , y = TMM ) ) +
        geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6 , aes(color = Type_new)) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 1 , aes(color = Type_new) ) +
        scale_color_manual(values=col) +
        facet_grid(.~Class) +
        xlab(NULL) +
        ylab("TMM")+
        theme_bw() +
        ggtitle(mut_type) +
        y_lab_lim +
        #geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black", face='bold') +
        stat_compare_means(comparisons = my_comparisons_1,method = "wilcox.test") +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='bottom',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold',hjust=0.5,vjust=0.5),
            strip.text.x = element_text(size = 12 , face = 'bold'),
            legend.text = element_text(size = 8,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.text.x = element_text(size = 10,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
            axis.ticks.length = unit(0.2, "cm") ,
            axis.line = element_line(size = 0.5))
    return(plot)
}

getTMM <- function(exp_use = exp_use , type = type){
    sample <- sapply( strsplit(names(exp_use) , "_") , "[" , 1)
    class <- sapply( strsplit(names(exp_use) , "_") , "[" , 2)
    class <- sapply( strsplit(class , "-") , "[" , 1)
    tmm <- as.numeric(exp_use)
    tmp_dat <- data.frame( Sample = sample , Class = class , TMM = tmm )
    ## 合并一个人的同一样本
    tmp_dat <- tmp_dat %>%
    group_by( Sample , Class ) %>%
    summarize( TMM = median(TMM) )
    tmp_dat <- data.frame(tmp_dat)
    tmp_dat$Type <- type

    return(tmp_dat)
}

describeGene <- function(dat_expression = dat_expression , image_path = image_path ){

    tmp_exp <- dat_expression

    ## 非共享突变型样本的表达
    type <- "Trunk"
    use_sample <- info_share
    exp_use <- tmp_exp[use_sample]
    tmm_mut_nonshare <- getTMM(exp_use = exp_use , type = type)

    ## 共享突变样本的表达
    type <- "Private"
    use_sample <- info_mut_nonshare
    exp_use <- tmp_exp[use_sample]
    tmm_mut_share <- getTMM(exp_use = exp_use , type = type)

    ## 野生型样本的表达
    type <- "Wild"
    use_sample <- info_wild
    exp_use <- tmp_exp[use_sample]
    tmm_wild <- getTMM(exp_use = exp_use , type = type)

    ## 合并
    tmm_combine <- rbind( tmm_wild , tmm_mut_nonshare , tmm_mut_share )
    print(tmm_combine$Class)
    tmm_combine$Class <- factor( tmm_combine$Class , levels = c( "All" , "IM" , "IGC" , "DGC" ) , order = T )
    tmm_combine$Type_new <- tmm_combine$Type

    #tmp_dat_use <- subset( tmm_combine , Class != "IM" & Type_new != "Wild" )
    tmp_dat_use <- subset( tmm_combine , Class != "IM" )
    if(length(unique(tmp_dat_use$Class)) > 1){
        from <- "NJSCC"
        y_tmm <- max(tmp_dat_use$TMM) * 1.1
        width <- 4.5
        title <- gene
        out_name <- paste0( image_path , "/" , gene , ".ShareVsMut." , mut_type , ".pdf" )
        plot <- plotTpm_combine(tmm_combine = tmp_dat_use , out_name = out_name , title = title , width = width , y_tmm = y_tmm)
        ggsave(file=out_name,plot=plot,width=width,height=4.5)

        tmp_dat_use <- tmp_dat_use %>%
        group_by( Type_new , Class ) %>%
        summarize( TMM = median(TMM) )
        out_name <- paste0( image_path , "/" , gene , ".ShareVsMut." , mut_type , ".tsv" )
        write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)
    }
}

##########################################################################################

dat_maf_public_use <- merge(dat_maf_public , info[,c("Tumor" , "Class_sub")] , by.x = "Sample" , by.y = "Tumor")
dat_maf_public_use$use_class <- paste0( dat_maf_public_use$ID , "_" , dat_maf_public_use$Class_sub  )
info$use_class <- paste0( info$ID , "_" , info$Class_sub )

## 突变型样本
mutTumor <- unique(dat_maf_public_use$use_class)
info_mut <- mutTumor
info_mut <- info_mut[info_mut %in% colnames(dat_expression)]

## 共享突变样本
shareSample <- subset( dat_maf_public_use , Share=="TRUE" )$use_class
info_share <- shareSample
info_share <- info_share[info_share %in% colnames(dat_expression)]

## 非共享突变样本
info_mut_nonshare <- info_mut[!(info_mut %in% info_share)]

##野生型样本
info_wild <- subset( info , !(use_class %in% mutTumor) )$use_class
info_wild <- info_wild[info_wild %in% colnames(dat_expression)]

## 分单个基因，看在一个人多个样本的改变情况
Hugo_Symbol <- gene
image_path <- out_path
dir.create(image_path , recursive = T)

describeGene(dat_expression = dat_expression , image_path = image_path )
